function ISStats = IS(X,Y,Tbeg,Tend,h,NWLags)
% lag/lead 
x = X(Tbeg:Tend-h); y = Y(Tbeg+h:Tend);
x = [ones(length(y),1) x];
[beta,tOls,~,tNwest,ISStats.R2,ISStats.adjr2] = ols(x,y,NWLags);
ISStats.beta = beta(2); ISStats.tOls = tOls(2); ISStats.tNwest = tNwest(2);

T = length(Y);
ISStats.ErrMdl = NaN(T,1); ISStats.ErrMdl(Tbeg+h:Tend) = y - x*beta;
ISStats.ErrUnc = NaN(T,1); ISStats.ErrUnc(Tbeg+h:Tend) = y - mean(y);

T = length(y);
Dum1 = zeros(T,1); Dum1(1:round(T/2)) = 1;
Dum2 = ones(T,1); Dum2(1:round(T/2)) = 0;
x = [Dum1 Dum2 x(:,2).*Dum1 x(:,2).*Dum2];
[ISStats.betaSplit,ISStats.tNwestSplit] = olssplit(x,y,NWLags);
return